logo

Όλα τα απαραίτητα αρχεία βρίσκονται στο github.

1 Πριν ξεκινήσουμε

Βεβαιώστε ότι το όνομα του project σας, όπως και το file path ΔΕΝ περιέχουν:
1. Ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Επί παραδείγματι το file path σας θα μπορούσε να είναι: C:/Environmental_Data_Analysis/Geographical_Data.

Όλα τα απαιτούμενα δεδομένα και αρχεία είναι ελεύθερα διαθέσιμα στο διαδίκτυο.

3 Διανυσματικά δεδομένα

3.1 Σύνορα χωρών

Θα χρησιμοποιήσουμε τα διανυσματικά δεδομένα της Ιταλίας, όπως έχουμε κάνει και σε προηγούμενα εργαστήρια.

Θυμηθείτε ότι θα χρειαστεί να αλλάξετε το file path και να χρησιμοποιήσετε τα δεδομένα από το Github.

4 Δεδομένα θέσης

4.1 Διαδικτυακή βάση δεδομένων GBIF

Θα χρησιμοποιήσουμε τα δεδομένα θέσης για όλα τα φυτικά είδη τα οποία απαντώνται στην Ιταλία, τα οποία περιέχονται στην διαδικτυακή βάση δεδομένων GBIF και από τα οποία έχουμε αφαιρέσει τις όποιες προβληματικές καταγραφές. Θα χρειαστεί να καθαρίσουμε τα ονόματα των taxa και στη συνέχεια να κρατήσουμε μόνο τα έγκυρα ονόματα.

## ============================================================================##
## Load the biotic cleaned data
## ============================================================================##
data_it_cleaned <- readRDS("Italian data cleaned.rds")
## ============================================================================##


## ============================================================================##
## First load the excel file with the correct species names
## ============================================================================##
nms <- readxl::read_excel("italian_names.xlsx")
## ============================================================================##


## ============================================================================##
## Keep only one instance per taxon from our occurrence data
## ============================================================================##
nms_flags <- data_it_cleaned %>% distinct(scientificName) %>% mutate(Taxon = nms$scientificName, 
    action = nms$action) %>% as_tibble()
## ============================================================================##


## ============================================================================##
## Then create the coords for the correct taxa
## ============================================================================##
coords <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action, 
    "DEL")) %>% dplyr::select(decimalLongitude, decimalLatitude, Taxon)

clean_dts <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action, 
    "DEL")) %>% distinct(Taxon, .keep_all = T)
## ============================================================================##

5 Αξιολόγηση κατά IUCN

Η αξιολόγηση του κινδύνου εξαφάνισης ενός είδους βασίζεται στην τυποποιημένη διαδικασία που έχει αναπτυχθεί από την Διεθνή Ένωση για την Προστασία της Φύσης (IUCN) και αποτελεί την πλέον αντικειμενική και αναλυτική προσέγγιση για την ανάδειξη αξόνων προτεραιότητας για την διατήρηση των ειδών και οικοτόπων. Ο Ερυθρός Κατάλογος των Απειλούμενων Ειδών ή IUCN Red List of Threatened Species παρέχει πληροφορίες σχετικά με την ταξινόμηση, την κατανομή, καθώς και το καθεστώς διατήρησης των φυτικών και ζωικών ειδών του πλανήτη μας, βάσει των κριτηρίων και των κατηγοριών της IUCN (IUCN Red List Categories and Criteria). Ο κύριος στόχος της διαδικασίας αυτής είναι να προσδιορίσει τον ενδεχόμενο κίνδυνο εξαφάνισης ενός είδους και κατά συνέπεια, να αναδείξει τα είδη εκείνα τα οποία χρήζουν άμεσης προτεραιότητας για την διατήρηση τους (Stévart et al., 2019).

Η δημιουργία ενός Ερυθρού Καταλόγου απαιτεί αξιόπιστα δεδομένα και προσεκτική εφαρμογή όλων των κριτηρίων της IUCN – μια διαδικασία η οποία είναι ιδιαίτερα χρονοβόρα. Το γεγονός αυτό αντανακλάται από το ότι υπάρχουν σημαντικά κενά στα έως τώρα αξιολογηθέντα είδη από την IUCN, καθώς μόλις 28.114 φυτικά είδη έχουν αξιολογηθεί μέχρι σήμερα. Το κενό αυτό είναι απόρροια της εξαιρετικά μεγάλης ποικιλότητας των φυτικών ειδών. Συνεπώς, παρότι η IUCN είναι κοντά στον στόχο που είχε θέσει, να έχει αξιολογήσει δηλαδή 38.500 φυτικά είδη έως το 2020 (Bachman et al., 2018), εντούτοις, ο αριθμός αυτός αποτελεί μόλις το 10% της φυτικής ποικιλότητας του πλανήτη και του δεύτερου στόχου της Παγκόσμιας Στρατηγικής για την Διατήρηση της Φυτικής ποικιλότητας (i.e., Global Strategy for Plant Conservation - περισσότερες πληροφορίες εδώ). Το γεγονός αυτό ώθησε την επιστημονική κοινότητα να προσπαθήσει να βρει τρόπους ώστε να επιταχύνει την διαδικασία αυτή.

6 Εκτίμηση κινδύνου εξαφάνισης - Εισαγωγικές έννοιες

Όπως είπαμε, μέχρι στιγμής έχει αξιολογηθεί ένα πολύ μικρό ποσοστό της έως τώρα γνωστής βιοποικιλότητας:

Η πλειονότητα των χωρών που αντιμετωπίζουν τα μεγαλύτερα προβλήματα όσον αφορά την κατάσταση διατήρησης των φυτικών taxa τα οποία φιλοξενούν εντοπίζονται στις λιγότερο αναπτυγμένες χώρες του πλανήτη:

Σε ευρωπαϊκό επίπεδο δε, γνωρίζουμε ότι η κατάσταση έχει προς το παρόν διαμορφωθεί ως εξής:

Όπου μπορούμε να διαπιστώσουμε ότι οι χώρες της Μεσογειακής λεκάνης και ιδιαίτερα η Ελλάδα, θεωρούνται θερμά σημεία κινδύνου εξαφάνισης:

Με την κυριότερη απειλή να είναι η εντατικοποίηση της βόσκησης

Όμως οι απειλές αυτές διαφέρουν τόσο μεταξύ των θερμών και των ψυχρών σημείων βιοποικιλότητας:

Όσο και στον χρόνο (για τα τελευταία 250 έτη τουλάχιστον):

Σύμφωνα με την IUCN, υπάρχουν 5 κατηγορίες κινδύνου εξαφάνισης2:

Ας τα δούμε πιο αναλυτικά:

Όσον αφορά δε την Ελλάδα:

Ενώ όσον αφορά την Κρήτη, το θερμότερο σημείο ενδημικής φυτικής ποικιλότητας στη Μεσόγειο:

Η κλιματική αλλαγή αναμένεται να έχει μεγάλες επιπτώσεις όσον αφορά τα θερμά σημεία ενδημικής φυτικής ποικιλότητας της νήσου (Kougioumoutzis et al., 2020a, 2020b)

7 Εκτίμηση κινδύνου εξαφάνισης - Στην πράξη

7.1 Φόρτωση των δεδομένων θέσης

Τώρα θα φορτώσουμε το set δεδομένων μας.

Στην περίπτωση που δεν έχετε κατεβάσει ήδη τα απαραίτητα αρχεία, κάντε το τώρα:

## [1] "Allium dilatatum"                    
## [2] "Campanula jacquinii subsp. jacquinii"
## [3] "Tulipa doerfleri"

7.2 Μορφοποιήση των δεδομένων

Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από την βιβλιοθήκη red. Θα χρειαστεί να φτιάξουμε τρια αντικείμενα, ένα για κάθε είδος το οποίο μας ενδιαφέρει και για το οποίο έχουμε δεδομένα. Σήμερα θα ασχοληθούμε τρια από τα τοπικά νησιωτικά ενδημικά της Κρήτης, το Allium dilalatum, την Campanula jacquinii subsp. jacquinii και την Tulipa doerfleri.

Το είδος Allium dilalatum - Photo from www.greekflora.gr

Το είδος Allium dilalatum - Photo from www.greekflora.gr

Το είδος Campanula jacquinii subsp. jacquinii - Photo from www.west-crete.com

Το είδος Campanula jacquinii subsp. jacquinii - Photo from www.west-crete.com

Το είδος Tulipa doerfleri- Photo from idr

Το είδος Tulipa doerfleri- Photo from idr

7.4 Εκτίμηση περιοχής κατάληψης (Area of Occupancy)

Τώρα είμαστε σε θέση να υπολογίσουμε την παρούσα περιοχή κατάληψης των ειδών που μας ενδιαφέρουν:

7.5 Εκτίμηση περιοχής εμφάνισης (Extent of Occurrence)

Τώρα είμαστε σε θέση να υπολογίσουμε την παρούσα περιοχή εμφάνισης των ειδών που μας ενδιαφέρουν:

7.7 Εκτίμηση περιοχής εμφάνισης (Extent of Occurrence) στο μέλλον

Ας υπολογίσουμε την παρούσα περιοχή εμφάνισης των ειδών που μας ενδιαφέρουν:

Τι παρατηρείτε; Υπάρχει διαφορά σε σχέση με την παρούσα κατάσταση; Πώς θα χαρακτηρίζατε τα είδη αυτά σύμφωνα με τις κατηγορίες κινδύνου εξαφάνισης της IUCN;

8 Ταχεία και μαζική αξιολόγηση ειδών κατά IUCN – The PACA approach

Η αξιολόγηση του καθεστώτος κινδύνου εξαφάνισης ενός είδους βάσει των κατευθυντήριων γραμμών της Διεθνούς Ενώσεως για την Διατήρηση της Φύσης IUCN αποτελούν τον πλέον αντικειμενικό και αναλυτικό τρόπο προσέγγισης για τον καθορισμό των στόχων προτεραιότητας διατήρησης. Ο ιστότοπος IUCN Red List of Threatened Species παρέχει πλειάδα πληροφοριών σχετικά με την ταξινομική, την κατανομή και το καθεστώς διατήρησης αρκετών φυτικών και ζωικών ειδών, βάσει των Κατηγοριών και των Κριτηρίων της IUCN. Ο κύριος στόχος αυτής της διαδικασίας είναι να προσδιοριστεί, χρησιμοποιώντας μια αναλυτικότατη και αντικειμενική μέθοδο, ο κίνδυνος εξαφάνισης ενός είδους και, συνεπώς, να προσδιοριστεί ποια είδη χρήζουν άμεσης προτεραιότητας διατήρησης. Ο κατάλογος Ερυθρών Δεδομένων της IUCN θεωρείται ως ο “χρυσός κανόνας” όσον αφορά τον εντοπισμό των ειδών εκείνων για τα οποία απαιτείται να δοθεί ιδιαίτερη προσοχή κατά τον σχεδιασμό και την υλοποίηση έργων περιβαλλοντικών επιπτώσεων (Stévart et al., 2019).

Η αξιολόγηση συνεπώς του καθεστώτος κινδύνου εξαφάνισης βάσει των Κριτηρίων της IUCN απαιτεί τόσο αξιόπιστα δεδομένα, όσο και προσεκτικότατη εφαρμογή και ενδελεχή γνώση των Κριτηρίων της IUCN και για αυτόν τον λόγο είναι μια ιδιαίτερα χρονοβόρα διαδικασία. Αυτό υποδηλώνεται και από το γεγονός ότι ένα μικρό ποσοστό των φυτικών ειδών (~10%) έχει έως τώρα αξιολογηθεί σύμφωνα με τα Κριτήρια της IUCN (Bachman et al., 2018) και από ότι φαίνεται ο Aichi Target 2 (Αξιολόγηση του καθεστώτος κινδύνου εξαφάνισης για την πλειονότητα των φυτικών ειδών - Global Strategy for Plant Conservation και εδώ), δεν πρόκειται να επιτευχθεί σύντομα. Για αυτόν τον λόγο, η επιστημονική κοινότητα καλέιται να επισπεύσει - στο μέτρο του δυνατού - την διαδικασία αξιολόγησης των φυτικών ειδών σε παγκόσμιο επίπεδο.

Πρόσφατα, προτάθηκε μια νέα προσέγγιση που βασίζεται στα κριτήρια της IUCN και παρέχει την δυνατότητα ταχείας, μαζικής και αξιόπιστης αξιολόγησης των ειδών μιας περιοχής/χώρας (Dauby et al., 2017; Stévart et al., 2019). Η προσέγγιση αυτή ονομάζεται Preliminary Automated Conservation Assessments (PACA) και αυτή θα χρησιμοποιήσουμε σήμερα.

Υπάρχει πιθανότητα αναλόγως των δυνατοτήτων του Η/Υ σας, η εντολή αυτή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή: italy_iucn_eval <- readRDS('Italian IUCN.rds').

Με αυτόν τον τρόπο λοιπόν, αξιολογήσαμε το καθεστώς κινδύνου εξαφάνισης όλων των ειδών που είχαμε κατεβάσει από την GBIF για την Ιταλία. Όμως, παρότι είχαμε προχωρήσει σε ορισμένα βήματα καθαρισμού των δεδομένων, απαιτούνται πολύ περισσότερα στάδια3 προκειμένου να αξιολογηθεί ένα είδος σύμφωνα με όλα τα κριτήρια της IUCN.

Μπορούμε επίσης να υπολογίσουμε και να αποθηκεύσουμε στον σκληρό μας δίσκο τα πολύγωνα της περιοχής εμφάνισης (Extent of Occurrence) για κάθε είδος, βασιζόμενοι σε δύο μεθόδους προσδιορισμού του πολυγώνου αυτού: την convex- και την alpha-hull μέθοδο (Εάν θέλετε, μπορείτε να διαβάσετε αυτήν την εργασία, εάν ενδιαφέρεστε να μάθετε κάτι περισσότερο σχετικά με τις μεθόδους αυτές).

Μην τρέξετε αυτή την εντολή τώρα.

8.1 Χωρικά πρότυπα του κινδύνου εξαφάνισης

Ας δούμε τώρα ποια περιοχή της Ιταλίας εμφανίζει το υψηλότερο ποσοστό κινδυνευόντων φυτικών ειδών:

##============================================================================##
## First create a common column and then join the occurrences with the
## evaluations
##============================================================================##
italy_iucn_eval <- italy_iucn_eval %>% rename(Taxon = taxa)
data_joined <- full_join(coords, italy_iucn_eval)
##============================================================================##


##============================================================================##
## Then select only the endangered taxa (CR, EN, VU)
##============================================================================##
endangered <- data_joined %>% filter(!str_detect(Category_code, 'LC'))
##============================================================================##


##============================================================================##
## Convert Italy to an sf spatial object
##============================================================================##
Italy_sf <- Italy %>% st_as_sf()
##============================================================================##


##============================================================================##
## Convert the endangered species to a sf spatial object
##============================================================================##
data_it_sf <- st_as_sf(endangered, 
                       coords = c("decimalLongitude", 
                                  "decimalLatitude"), 
                       crs = 4326)
##============================================================================##


##============================================================================##
## Make sure that all these occurrences lie within Greece
##============================================================================##
point_endangered <- data_it_sf[Italy_sf,]
##============================================================================##


##============================================================================##
## Intersect Greece with these occurrences
##============================================================================##
endangered_intersection <- Italy_sf %>%
  st_join(point_endangered) %>%
  group_by(NOME_PRO) %>% ## Here we group the data based on the county names
  
  ## First we create a variable containing the sum of the species in each county
  summarize(Species_Number = sum(NROW(unique(Taxon))),
            
            ## Then we just sum the number of occurrences in each county
            Occurrence_Number = n(),
            
            ## Finally, we calculate the proportion of the endangered species in
            ## each county
            Risk = sum(NROW(unique(Taxon)))/sum(NROW(unique(data_joined$Taxon)))*100
  )
##============================================================================##


##============================================================================##
## Colour function
##============================================================================##
col2alpha <- function(col, alpha) {
  col_rgb <- col2rgb(col)/255
  rgb(col_rgb[1], col_rgb[2], col_rgb[3], 
      alpha = alpha)
}
##============================================================================##
##============================================================================##
## Create the main map
##============================================================================##
main_map <- ggplot(Italy_sf) + 
  geom_sf(fill = "antiquewhite1") + 
  geom_sf() + 
  geom_sf(data = Italy_sf, 
          fill = NA) + 
  geom_sf(aes(fill = Risk), 
          data = endangered_intersection, 
          inherit.aes = FALSE) + 
  theme(axis.title = element_blank(),
        legend.position = "bottom") +
  coord_sf(expand = FALSE) +
  scale_fill_viridis_c(option = "viridis",
                       trans = "sqrt",
                       limits = c(0, 12.4),
                       breaks = c(0.4, 3, 12.3),
                       labels = scales::label_percent(scale = 1),
                       name = expression(bold('Extinction risk'))) +
  geom_rect(
    xmin = st_bbox(Italy_sf_tr)[[1]],
    ymin = st_bbox(Italy_sf_tr)[[2]],
    xmax = st_bbox(Italy_sf_tr)[[3]],
    ymax = st_bbox(Italy_sf_tr)[[4]],
    fill = NA, 
    colour = "red",
    size = 0.6
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.25) +
  annotation_north_arrow(location = "bl", 
                         which_north = "true", 
                         pad_x = unit(0.75, "in"), 
                         pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  geom_sf(data = Italian_counties,
          fill = NA,
          lwd = 1,
          color = 'black') + 
  geom_sf_label_repel(aes(label = NOME_REG),
                     data = Italian_counties) + 
  theme(panel.grid.major = element_line(color = gray(.5),
                                        linetype = "dashed", 
                                        size = 0.5), 
        panel.background = element_rect(fill = col2alpha('steelblue', 0.2))) +
  labs(title = 'Extinction risk in Italian counties',
       subtitle = 'Source: GBIF')
##============================================================================##


##============================================================================##
## Create the inset map
##============================================================================##
side_map <- ggplot(Italy_sf) + 
  geom_sf(fill = "antiquewhite1") + 
  geom_sf() + 
  geom_sf(data = Italy_sf, 
          fill = NA) + 
  geom_sf(aes(fill = Risk), 
          data = endangered_intersection, 
          inherit.aes = FALSE) + 
  theme(axis.title = element_blank(),
        legend.position = "bottom") +
  coord_sf(expand = FALSE) +
  scale_fill_viridis_c(option = "viridis",
                       trans = "sqrt",
                       limits = c(0, 12.4),
                       breaks = c(0.4, 3, 12.3),
                       labels = scales::label_percent(scale = 1),
                       name = 'Extinction risk') +
  geom_rect(
    xmin = st_bbox(Italy_sf_tr)[[1]],
    ymin = st_bbox(Italy_sf_tr)[[2]],
    xmax = st_bbox(Italy_sf_tr)[[3]],
    ymax = st_bbox(Italy_sf_tr)[[4]],
    fill = NA, 
    colour = "red",
    size = 0.6
  ) +
  geom_sf(data = Italian_counties,
          fill = NA,
          lwd = 1,
          color = 'black') + 
  geom_sf_label_repel(aes(label = NOME_PRO),
                      data = Italy_sf_tr) +
  coord_sf(xlim = c(st_bbox(Italy_sf_tr)[[1]],
                    st_bbox(Italy_sf_tr)[[3]]),
           ylim = c(st_bbox(Italy_sf_tr)[[2]],
                    st_bbox(Italy_sf_tr)[[4]]),
           expand = FALSE)  + 
  theme_void() +
  theme(legend.position = "none")
##============================================================================##


##============================================================================##
## And the final map
##============================================================================##
main_map %>% 
  ggdraw() +
  draw_plot(side_map,
            
    # The distance along a (0,1) x-axis to draw the left edge of the plot
    x = 0.625, 
    
    # The distance along a (0,1) y-axis to draw the bottom edge of the plot
    y = 0.65,
    
    # The width and height of the plot expressed as proportion of the entire
    # ggdraw object
    width = 0.35, 
    height = 0.35)

Ποια περιοχή της Ιταλίας εμφανίζει το υψηλότερο ποσοστό κινδυνευόντων ειδών;

NOME_PRO Species_Number Occurrence_Number Risk geometry
BOLZANO/BOZEN 783 2108 12.3990499 MULTIPOLYGON (((12.20557 47…
COMO 467 2585 7.3950911 MULTIPOLYGON (((9.279479 46…
VALLE D’AOSTA 400 1020 6.3341251 MULTIPOLYGON (((7.58858 45….
VERBANO-CUSIO-OSSOLA 392 1409 6.2074426 MULTIPOLYGON (((8.449767 46…
SONDRIO 347 1118 5.4948535 MULTIPOLYGON (((10.2493 46….
L’AQUILA 345 695 5.4631829 MULTIPOLYGON (((13.40447 42…
TORINO 302 467 4.7822644 MULTIPOLYGON (((7.859052 45…
TRENTO 279 413 4.4180523 MULTIPOLYGON (((11.82411 46…
CUNEO 267 495 4.2280285 MULTIPOLYGON (((7.990904 44…
VARESE 266 884 4.2121932 MULTIPOLYGON (((8.78101 46….
PALERMO 264 403 4.1805226 MULTIPOLYGON (((13.07012 38…
GROSSETO 214 639 3.3887569 MULTIPOLYGON (((11.13574 42…
BERGAMO 199 301 3.1512272 MULTIPOLYGON (((10.10323 46…
BELLUNO 184 260 2.9136975 MULTIPOLYGON (((12.50593 46…
FOGGIA 169 392 2.6761679 MULTIPOLYGON (((16.06218 41…
UDINE 166 250 2.6286619 MULTIPOLYGON (((12.76338 46…
NUORO 157 283 2.4861441 MULTIPOLYGON (((9.583129 40…
BRESCIA 150 188 2.3752969 MULTIPOLYGON (((10.49807 46…
SIRACUSA 150 293 2.3752969 MULTIPOLYGON (((15.11665 37…
ALESSANDRIA 138 222 2.1852732 MULTIPOLYGON (((8.404617 45…
CATANIA 138 270 2.1852732 MULTIPOLYGON (((15.09013 37…
OGLIASTRA 136 198 2.1536025 MULTIPOLYGON (((9.628694 40…
TRAPANI 128 175 2.0269200 MULTIPOLYGON (((12.44316 37…
REGGIO DI CALABRIA 112 136 1.7735550 MULTIPOLYGON (((15.72781 37…
MESSINA 109 135 1.7260491 MULTIPOLYGON (((14.96671 38…
SIENA 104 154 1.6468725 MULTIPOLYGON (((11.42858 43…
ROMA 101 160 1.5993666 MULTIPOLYGON (((12.62241 41…
FROSINONE 97 178 1.5360253 MULTIPOLYGON (((13.33753 41…
BOLOGNA 95 144 1.5043547 MULTIPOLYGON (((11.29036 44…
TREVISO 91 106 1.4410135 MULTIPOLYGON (((12.34077 46…
GENOVA 90 125 1.4251781 MULTIPOLYGON (((8.686318 44…
VERONA 84 110 1.3301663 MULTIPOLYGON (((10.88381 45…
VERCELLI 83 158 1.3143310 MULTIPOLYGON (((8.204472 45…
PERUGIA 82 149 1.2984956 MULTIPOLYGON (((12.36166 43…
SASSARI 80 133 1.2668250 MULTIPOLYGON (((8.203562 40…
LIVORNO 78 105 1.2351544 MULTIPOLYGON (((10.54294 42…
RIETI 76 102 1.2034838 MULTIPOLYGON (((13.28007 42…
IMPERIA 75 114 1.1876485 MULTIPOLYGON (((7.760733 43…
OLBIA-TEMPIO 75 112 1.1876485 MULTIPOLYGON (((9.515061 40…
TRIESTE 75 101 1.1876485 MULTIPOLYGON (((13.69162 45…
CAGLIARI 70 81 1.1084719 MULTIPOLYGON (((9.309161 39…
MILANO 67 99 1.0609660 MULTIPOLYGON (((8.93997 45….
PESCARA 65 78 1.0292953 MULTIPOLYGON (((14.07289 42…
CHIETI 64 89 1.0134600 MULTIPOLYGON (((14.30305 42…
FIRENZE 64 77 1.0134600 MULTIPOLYGON (((11.44765 44…
CALTANISSETTA 62 98 0.9817894 MULTIPOLYGON (((13.81375 37…
CARBONIA-IGLESIAS 62 104 0.9817894 MULTIPOLYGON (((8.361779 39…
SAVONA 62 79 0.9817894 MULTIPOLYGON (((8.192667 44…
LA SPEZIA 61 82 0.9659541 MULTIPOLYGON (((9.567304 44…
LECCO 60 97 0.9501188 MULTIPOLYGON (((9.415008 46…
TARANTO 60 122 0.9501188 MULTIPOLYGON (((17.29291 40…
POTENZA 59 85 0.9342835 MULTIPOLYGON (((15.72353 39…
LECCE 57 99 0.9026128 MULTIPOLYGON (((18.27608 39…
LUCCA 55 74 0.8709422 MULTIPOLYGON (((10.16621 43…
AGRIGENTO 48 89 0.7600950 MULTIPOLYGON (((13.02734 37…
ORISTANO 48 68 0.7600950 MULTIPOLYGON (((8.497556 40…
TERAMO 48 84 0.7600950 MULTIPOLYGON (((13.92029 42…
BRINDISI 47 59 0.7442597 MULTIPOLYGON (((17.96702 40…
MATERA 47 65 0.7442597 MULTIPOLYGON (((16.2523 40….
PADOVA 47 54 0.7442597 MULTIPOLYGON (((11.83986 45…
VICENZA 47 87 0.7442597 MULTIPOLYGON (((11.53894 46…
PARMA 43 53 0.6809184 MULTIPOLYGON (((10.1873 45….
COSENZA 42 50 0.6650831 MULTIPOLYGON (((15.81585 39…
RAGUSA 42 73 0.6650831 MULTIPOLYGON (((14.49321 36…
VENEZIA 42 43 0.6650831 MULTIPOLYGON (((12.79997 45…
AVELLINO 40 49 0.6334125 MULTIPOLYGON (((15.19857 41…
BARI 40 57 0.6334125 MULTIPOLYGON (((16.65049 41…
ISERNIA 39 83 0.6175772 MULTIPOLYGON (((14.30542 41…
ASTI 38 40 0.6017419 MULTIPOLYGON (((8.046812 45…
SALERNO 38 42 0.6017419 MULTIPOLYGON (((15.32929 40…
FORLI’-CESENA 37 44 0.5859066 MULTIPOLYGON (((12.07366 44…
MACERATA 37 54 0.5859066 MULTIPOLYGON (((13.72925 43…
MEDIO CAMPIDANO 37 44 0.5859066 MULTIPOLYGON (((8.950374 39…
NAPOLI 35 40 0.5542359 MULTIPOLYGON (((14.03959 40…
RAVENNA 35 47 0.5542359 MULTIPOLYGON (((12.28071 44…
PAVIA 34 40 0.5384006 MULTIPOLYGON (((8.847651 45…
VITERBO 34 56 0.5384006 MULTIPOLYGON (((11.86313 42…
PISA 32 42 0.5067300 MULTIPOLYGON (((10.42773 43…
PORDENONE 29 40 0.4592241 MULTIPOLYGON (((12.52112 46…
LATINA 28 30 0.4433888 MULTIPOLYGON (((12.94157 41…
CATANZARO 26 28 0.4117181 MULTIPOLYGON (((16.62624 39…
PESARO E URBINO 26 27 0.4117181 MULTIPOLYGON (((12.76836 43…
ROVIGO 26 35 0.4117181 MULTIPOLYGON (((12.32842 45…
FERRARA 25 39 0.3958828 MULTIPOLYGON (((12.01548 44…
MODENA 25 27 0.3958828 MULTIPOLYGON (((11.08931 44…
AREZZO 24 33 0.3800475 MULTIPOLYGON (((11.72191 43…
CREMONA 24 26 0.3800475 MULTIPOLYGON (((9.53991 45….
ENNA 24 63 0.3800475 MULTIPOLYGON (((14.57105 37…
MANTOVA 24 27 0.3800475 MULTIPOLYGON (((10.66871 45…
ANCONA 23 30 0.3642122 MULTIPOLYGON (((13.50006 43…
ASCOLI PICENO 21 25 0.3325416 MULTIPOLYGON (((13.78665 43…
GORIZIA 20 21 0.3167063 MULTIPOLYGON (((13.53041 45…
VIBO VALENTIA 20 22 0.3167063 MULTIPOLYGON (((15.90041 38…
MONZA E DELLA BRIANZA 19 28 0.3008709 MULTIPOLYGON (((9.253944 45…
MASSA-CARRARA 17 21 0.2692003 MULTIPOLYGON (((9.903351 44…
PIACENZA 16 21 0.2533650 MULTIPOLYGON (((9.917413 45…
CASERTA 15 16 0.2375297 MULTIPOLYGON (((14.1601 41….
TERNI 15 15 0.2375297 MULTIPOLYGON (((12.1291 42….
BIELLA 14 15 0.2216944 MULTIPOLYGON (((8.129534 45…
NOVARA 13 14 0.2058591 MULTIPOLYGON (((8.496885 45…
PISTOIA 13 14 0.2058591 MULTIPOLYGON (((10.64794 44…
LODI 11 12 0.1741884 MULTIPOLYGON (((9.442017 45…
PRATO 11 12 0.1741884 MULTIPOLYGON (((11.1729 44….
CAMPOBASSO 8 8 0.1266825 MULTIPOLYGON (((14.84675 42…
REGGIO NELL’EMILIA 7 8 0.1108472 MULTIPOLYGON (((10.72925 44…
BARLETTA-ANDRIA-TRANI 5 5 0.0791766 MULTIPOLYGON (((15.95552 41…
BENEVENTO 4 4 0.0633413 MULTIPOLYGON (((15.0303 41….
CROTONE 4 11 0.0633413 MULTIPOLYGON (((17.05156 39…
FERMO 4 4 0.0633413 MULTIPOLYGON (((13.74288 43…
RIMINI 4 4 0.0633413 MULTIPOLYGON (((12.49431 44…

9 Evolutionary Distinct and Globally Endangered

Όπως έχουμε αναφέρει σε προηγούμενο εργαστήριο, στο φυλογενετικό επίπεδο, οι ισοδύναμοι των παραδοσιακών, ταξινομικών δεικτών για τον εντοπισμό θερμών σημείων βιοποικιλότητας (ΘΣΒ), είναι η φυλογενετική ποικιλότητα [ΦΠ, sensu Faith (1992) – ισοδύναμη του αριθμού των ειδών], ο φυλογενετικός ενδημισμός [ΦΕ, sensu Rosauer et al. (2009) – ισοδύναμος του αριθμού των ενδημικών/ποσοστού ενδημισμού] και ο δείκτης EDGE [Evolutionary Distinct and Globally Endangered – EDGE, sensu Isaac et al. (2007) – ισοδύναμος της πιθανότητας εξαφάνισης]. Οι εν λόγω δείκτες ποσοτικοποιούν διαφορετικές πτυχές της εξελικτικής ποικιλότητας (Tucker et al., 2017).

Η ΦΠ είναι το άθροισμα του μήκους των κλάδων που συνδέουν ένα σύμπλεγμα ειδών με την ρίζα ενός φυλογενετικού δένδρου (Faith, 1992).

Ο ΦΕ ποσοτικοποιεί τον βαθμό στον οποίο ένα σημαντικό ποσοστό της ΦΠ εντοπίζεται αποκλειστικά εντός της εκάστοτε περιοχής μελέτης (Rosauer et al., 2009).

Ο δείκτης EDGE συνδυάζει την εξελικτική διακριτότητα (ED, η φυλογενετική απομόνωση ενός είδους) με το καθεστώς κινδύνου εξαφάνισης σε παγκόσμιο επίπεδο (GE) σύμφωνα με τα πρότυπα της IUCN (Isaac et al., 2007).

Οι δείκτες ΦΕ και EDGE αντιπροσωπεύουν σταθμισμένες παραλλαγές της ΦΠ σε γεωγραφική κλίμακα και καθεστώς απειλής, αντίστοιχα (Daru et al., 2019).

Συνεπώς ο δείκτης EDGE μπορεί να αποβεί ιδιαίτερα χρήσιμος κατά τον εντοπισμό και την προτεραιοποίηση των θερμών σημείων βιοποικιλότητας, ανεξαρτήτως κλίμακας (Daru et al., 2019).

Σήμερα λοιπόν, θα υπολογίσουμε τον δείκτη EDGE για τα φυτικά taxa τα οποία απαντώνται στην Ιταλία, σύμφωνα με τα δεδομένα τα οποία είναι διαθέσιμα στην βάση δεδομένων GBIF. θα πρέπει να τονιστεί ότι ο δείκτης EDGE πρέπει να υπολογίζεται για τα ενδημικά taxa της περιοχής μελέτης και όχι για το σύνολο των ειδών που απαντώνται σε αυτήν.

9.1 Υπολογισμός του δείκτη EDGE

Ας υπολογίσουμε τον δείκτη EDGE για τα είδη τα οποία απαντώνται στην Ιταλία.

## [1] "Cycadeoidea_etrusca"

Εάν σας δείτε κάποιο error στην επόμενη εντολή, τρέξτε τις σειρές του κώδικα που εμφανίζονται ως σχόλια.

Ας δούμε ποια είναι η κατανομή των τιμών του δείκτη EDGE για τα φυτικά taxa της Ιταλίας, κάποια στατιστικά στοιχεία για τον δείκτη αυτόν, καθώς και ποια είδη (αλλά και γένη και οικογένειες) εμφανίζουν τις υψηλότερες τιμές για τον δείκτη EDGE και κατά συνέπεια θα πρέπει να δώσουμε μεγαλύτερη έμφαση στην προτεραιοποίηση τους:

Data summary
Name italian_edge_tbl
Number of rows 6314
Number of columns 5
_______________________
Column type frequency:
character 4
numeric 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
order 0 1 6 15 0 59 0
family 0 1 7 16 0 179 0
genus 0 1 3 20 0 1222 0
Taxon 0 1 8 48 0 6314 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
EDGE 0 1 3.02 0.86 0.64 2.4 2.94 3.57 6.32 ▁▇▇▂▁
family genus Mean_EDGE Species_Number
Orchidaceae Ophrys 3.052281 122
Cyperaceae Carex 2.166766 97
Asteraceae Centaurea 2.121417 93
Asteraceae Hieracium 2.114490 84
Euphorbiaceae Euphorbia 2.933861 70
Caryophyllaceae Silene 2.988936 66
Taxon EDGE
Huperzia_selago 6.315989
Spinulum_annotinum 6.087944
Equisetum_hyemale_subsp_hyemale 6.055210
Equisetum_meridionale 6.055210
Cycas_revoluta 6.036176
Equisetum_arvense 5.789793

9.2 Χωρικά πρότυπα του δείκτη EDGE - Θερμά σημεία βιοποικιλότητας

Τώρα μπορούμε να υπολογίσουμε ποιες περιοχές της Ιταλίας εμφανίζουν υψηλές τιμές για τον δείκτη EDGE και μπορούν να θεωρηθούν ως θερμά ή ψυχρά σημεία βιοποικιλότητας. Πρώτα θα χρειαστεί να εισάγαγουμε το αρχείο το οποίο είχαμε δημιουργήσει σε προηγούμενο εργαστήριο και αφορά τη μήτρα παρουσίας/απουσίας των φυτικών taxa στην Ιταλία:

Και στη συνέχεια, μπορούμε να δημιουργήσουμε ένα χωρικό αντικείμενο το οποίο θα περιέχει πληροφορία σχετικά με το ποιες περιοχές δρουν ως θερμά ή ψυχρά σημεία βιοποικιλότητας (ταξινομικής και φυλογενετικής), καθώς και ως θερμά ή ψυχρά σημεία κινδυνού εξαφάνισης (δείκτης EDGE):

10 Appendix: R code

##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
options(Encoding="UTF-8")

library(knitr)
library(rmdformats)

## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
                 cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               # comment=NA,
               message=FALSE,
               warning=FALSE,
               eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
htmltools::img(src = knitr::image_uri(file.path("G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/HBS Seminar/Logo_HBS/ebe-LOGO-005-300_1_70.jpg")), 
               alt = 'logo', 
               style = 'position:absolute; top:0; right:0; padding:5px;')
##===========================================================================##
## Load them
##===========================================================================##
library(raster)
library(sf)
library(tidyverse)
library(ConR)
library(red)
library(rasterVis)
library(ggspatial)
library(picante)
library(caper)
library(phyloregion)
library(ggsflabel)
library(cowplot)
##===========================================================================##
##===========================================================================##
## Load the spatial data for Italian provinces
##===========================================================================##
Italy <- read_sf('Italin provinces.shp') %>% 
  st_transform(4326) %>% 
  as_Spatial()
##===========================================================================##


##===========================================================================##
## Load the spatial data for Italian counties
##===========================================================================##

Italian_counties <- read_sf("reg2011_g.shp") %>% 
  mutate(NOME_REG = str_replace_all(NOME_REG, 
                                    c("VALLE D'AOSTA/VALLÉE D'AOSTE\r\nVALLE D'AOSTA/VALLÉE D'AOSTE" = "AOSTA",
                                                "TRENTINO-ALTO ADIGE/SUDTIROL" = "TRENTINO",
                                                "FRIULI VENEZIA GIULIA" ="FRIULI",
                                                "EMILIA-ROMAGNA" ="EMILIA"))) %>% 
  st_transform(4326)
##===========================================================================##


##===========================================================================##
## Keep only the names in Trentino
##===========================================================================##
Italy_sf_tr <- Italy %>% 
  st_as_sf() %>% 
  st_join(Italian_counties) %>% 
  filter(str_detect(NOME_REG, 'TREN'))
##===========================================================================##
##============================================================================##
## Load the biotic cleaned data                             
##============================================================================##
data_it_cleaned <- readRDS('Italian data cleaned.rds')
##============================================================================##


##============================================================================##
## First load the excel file with the correct species names
##============================================================================##
nms <- readxl::read_excel('italian_names.xlsx')
##============================================================================##


##============================================================================##
## Keep only one instance per taxon from our occurrence data
##============================================================================##
nms_flags <- data_it_cleaned %>% 
  distinct(scientificName) %>% 
  mutate(Taxon = nms$scientificName,
         action = nms$action) %>% 
  as_tibble() 
##============================================================================##


##============================================================================##
## Then create the coords for the correct taxa
##============================================================================##
coords <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(decimalLongitude, 
                decimalLatitude,
                Taxon)

clean_dts <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  distinct(Taxon, .keep_all = T)
##============================================================================##

pre {
  max-height: 400px;
  overflow-y: auto;
}

pre[class] {
  max-height: 400px;
}
library(magick)
nms <- list.files('MAGICK',
                  full.names = T)

lapply(nms, image_read) %>% 
  image_join() %>% 
  image_animate(fps = 1)
aegean %>% 
  download_this(
    output_name = "PVA",
    output_extension = ".xlsx",
    button_label = "Download the dataset",
    button_type = "success"
  )
library(downloadthis)
aegean <- readxl::read_excel("./PVA_Crete.xlsx")
##==================================================================================##
## Load the data
##==================================================================================##
cretan_taxa <- readxl::read_excel("./PVA_Crete.xlsx") %>% 
  dplyr::rename(decimalLongitude = lon,
                decimalLatitude = lat)

cretan_taxa$Taxon <- factor(cretan_taxa$Taxon)
levels(cretan_taxa$Taxon)

##==================================================================================##
## Create 3 objects
##==================================================================================##
all_dila <- cretan_taxa %>% filter(Taxon == 'Allium dilatatum')
camp_jac <- cretan_taxa %>% filter(Taxon == 'Campanula jacquinii subsp. jacquinii')
tul_doer <- cretan_taxa %>% filter(Taxon == 'Tulipa doerfleri')
##==================================================================================##
##==================================================================================##
## Convert to matrices
##==================================================================================##
all_dila_mat <- all_dila %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
camp_jac_mat <- camp_jac %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
tul_doer_mat <- tul_doer %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
##==================================================================================##
##==================================================================================##
## Area of occupancy
##==================================================================================##
aoo(all_dila_mat) 
aoo(camp_jac_mat) 
aoo(tul_doer_mat) 
##==================================================================================##
##==================================================================================##
## Extent of Occurrence
##==================================================================================##
eoo(all_dila_mat) 
eoo(camp_jac_mat) 
eoo(tul_doer_mat) 
##==================================================================================##
##==================================================================================##
## Load the rds files
##==================================================================================##
allium_future <- readRDS('./RDS/Allium future.rds')
campanula_future <- readRDS('./RDS/Campanula future.rds')
tulipa_future <- readRDS('./RDS/Tulipa future.rds')
##==================================================================================##


##==================================================================================##
## Area of occupancy
##==================================================================================##
aoo(allium_future) 
aoo(campanula_future) 
aoo(tulipa_future) 
##==================================================================================##

##==================================================================================##
## Extent of Occurrence
##==================================================================================##
eoo(allium_future) 
eoo(campanula_future) 
eoo(tulipa_future) 
##==================================================================================##
##============================================================================##
## First, let's select only the required columns
##============================================================================##
data_it_sel <- coords %>% 
  dplyr::select(decimalLatitude,
                decimalLongitude,
                Taxon)
##============================================================================##


##============================================================================##
## Then, let's evaluate our taxa
##============================================================================##
italy_iucn_eval <- IUCN.eval(data_it_sel,  
                                
                              ## Here we use the spatial object for the country
                              ## we want
                              country_map = Italy,  
                                
                              ## We can run the function in parallel to save
                              ## time
                              parallel = T,  
                                
                              ## Change the number of cores accordingly
                              NbeCores = 23 ) %>%   
    
  ## Here we create a new column (EOO_tr) and we substitute the NA values with
  ## the values from the AOO column (Why?)
  mutate(EOO_tr = ifelse(is.na(EOO), 
                       AOO, 
                       EOO)) %>% 
  as_tibble() %>% 
  dplyr::select(taxa, 
                EOO_tr, 
                everything())


##============================================================================##
## Let's see the result
##============================================================================##
italy_iucn_eval
##============================================================================##
## First we run the convex-hull method
##============================================================================##
eooshp_convexhull <- EOO.computing(data_gr_sel, 
                        exclude.area = T,
                        country_map = Greece,
                        export_shp = T,
                        write_shp = T,
                        Name_Sp = 'scientificName',
                        parallel = T,
                        NbeCores = 4, ## Change them
                        show_progress = T,
                        method.less.than3 = 'arbitrary') ## What does this option mean?

convexhull <- eooshp_convexhull$spatial.polygon_1
##============================================================================##


##============================================================================##
## Then the alpha-hull method
##============================================================================##
eooshp_ahull <- EOO.computing(data_gr_sel, 
                         exclude.area = T,
                         country_map = Greece,
                         export_shp = T,
                         write_shp = T,
                         Name_Sp = 'scientificName',
                         parallel = T,
                         NbeCores = 4, ## Change them
                         show_progress = T,
                         method.less.than3 = 'arbitrary',
                         method.range = 'alpha.hull')

ahull <- eooshp_ahull$spatial.polygon_1
##============================================================================##
##============================================================================##
## First create a common column and then join the occurrences with the
## evaluations
##============================================================================##
italy_iucn_eval <- italy_iucn_eval %>% rename(Taxon = taxa)
data_joined <- full_join(coords, italy_iucn_eval)
##============================================================================##


##============================================================================##
## Then select only the endangered taxa (CR, EN, VU)
##============================================================================##
endangered <- data_joined %>% filter(!str_detect(Category_code, 'LC'))
##============================================================================##


##============================================================================##
## Convert Italy to an sf spatial object
##============================================================================##
Italy_sf <- Italy %>% st_as_sf()
##============================================================================##


##============================================================================##
## Convert the endangered species to a sf spatial object
##============================================================================##
data_it_sf <- st_as_sf(endangered, 
                       coords = c("decimalLongitude", 
                                  "decimalLatitude"), 
                       crs = 4326)
##============================================================================##


##============================================================================##
## Make sure that all these occurrences lie within Greece
##============================================================================##
point_endangered <- data_it_sf[Italy_sf,]
##============================================================================##


##============================================================================##
## Intersect Greece with these occurrences
##============================================================================##
endangered_intersection <- Italy_sf %>%
  st_join(point_endangered) %>%
  group_by(NOME_PRO) %>% ## Here we group the data based on the county names
  
  ## First we create a variable containing the sum of the species in each county
  summarize(Species_Number = sum(NROW(unique(Taxon))),
            
            ## Then we just sum the number of occurrences in each county
            Occurrence_Number = n(),
            
            ## Finally, we calculate the proportion of the endangered species in
            ## each county
            Risk = sum(NROW(unique(Taxon)))/sum(NROW(unique(data_joined$Taxon)))*100
  )
##============================================================================##


##============================================================================##
## Colour function
##============================================================================##
col2alpha <- function(col, alpha) {
  col_rgb <- col2rgb(col)/255
  rgb(col_rgb[1], col_rgb[2], col_rgb[3], 
      alpha = alpha)
}
##============================================================================##
##============================================================================##
## Create the main map
##============================================================================##
main_map <- ggplot(Italy_sf) + 
  geom_sf(fill = "antiquewhite1") + 
  geom_sf() + 
  geom_sf(data = Italy_sf, 
          fill = NA) + 
  geom_sf(aes(fill = Risk), 
          data = endangered_intersection, 
          inherit.aes = FALSE) + 
  theme(axis.title = element_blank(),
        legend.position = "bottom") +
  coord_sf(expand = FALSE) +
  scale_fill_viridis_c(option = "viridis",
                       trans = "sqrt",
                       limits = c(0, 12.4),
                       breaks = c(0.4, 3, 12.3),
                       labels = scales::label_percent(scale = 1),
                       name = expression(bold('Extinction risk'))) +
  geom_rect(
    xmin = st_bbox(Italy_sf_tr)[[1]],
    ymin = st_bbox(Italy_sf_tr)[[2]],
    xmax = st_bbox(Italy_sf_tr)[[3]],
    ymax = st_bbox(Italy_sf_tr)[[4]],
    fill = NA, 
    colour = "red",
    size = 0.6
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.25) +
  annotation_north_arrow(location = "bl", 
                         which_north = "true", 
                         pad_x = unit(0.75, "in"), 
                         pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  geom_sf(data = Italian_counties,
          fill = NA,
          lwd = 1,
          color = 'black') + 
  geom_sf_label_repel(aes(label = NOME_REG),
                     data = Italian_counties) + 
  theme(panel.grid.major = element_line(color = gray(.5),
                                        linetype = "dashed", 
                                        size = 0.5), 
        panel.background = element_rect(fill = col2alpha('steelblue', 0.2))) +
  labs(title = 'Extinction risk in Italian counties',
       subtitle = 'Source: GBIF')
##============================================================================##


##============================================================================##
## Create the inset map
##============================================================================##
side_map <- ggplot(Italy_sf) + 
  geom_sf(fill = "antiquewhite1") + 
  geom_sf() + 
  geom_sf(data = Italy_sf, 
          fill = NA) + 
  geom_sf(aes(fill = Risk), 
          data = endangered_intersection, 
          inherit.aes = FALSE) + 
  theme(axis.title = element_blank(),
        legend.position = "bottom") +
  coord_sf(expand = FALSE) +
  scale_fill_viridis_c(option = "viridis",
                       trans = "sqrt",
                       limits = c(0, 12.4),
                       breaks = c(0.4, 3, 12.3),
                       labels = scales::label_percent(scale = 1),
                       name = 'Extinction risk') +
  geom_rect(
    xmin = st_bbox(Italy_sf_tr)[[1]],
    ymin = st_bbox(Italy_sf_tr)[[2]],
    xmax = st_bbox(Italy_sf_tr)[[3]],
    ymax = st_bbox(Italy_sf_tr)[[4]],
    fill = NA, 
    colour = "red",
    size = 0.6
  ) +
  geom_sf(data = Italian_counties,
          fill = NA,
          lwd = 1,
          color = 'black') + 
  geom_sf_label_repel(aes(label = NOME_PRO),
                      data = Italy_sf_tr) +
  coord_sf(xlim = c(st_bbox(Italy_sf_tr)[[1]],
                    st_bbox(Italy_sf_tr)[[3]]),
           ylim = c(st_bbox(Italy_sf_tr)[[2]],
                    st_bbox(Italy_sf_tr)[[4]]),
           expand = FALSE)  + 
  theme_void() +
  theme(legend.position = "none")
##============================================================================##


##============================================================================##
## And the final map
##============================================================================##
main_map %>% 
  ggdraw() +
  draw_plot(side_map,
            
    # The distance along a (0,1) x-axis to draw the left edge of the plot
    x = 0.625, 
    
    # The distance along a (0,1) y-axis to draw the bottom edge of the plot
    y = 0.65,
    
    # The width and height of the plot expressed as proportion of the entire
    # ggdraw object
    width = 0.35, 
    height = 0.35)

##============================================================================##
## The solution
##============================================================================##
endangered_intersection %>% arrange(desc(Risk))
##============================================================================##
## Load the tree
##============================================================================##
italy_tree <- readRDS('Italian phylogenetic tree.rds')


##============================================================================##
## See which taxa are missing
##============================================================================##
setdiff(italy_iucn_eval$Taxon, italy_tree$tip.label)
##============================================================================##
## Create the columns needed for the estimation of the EDGE index
##============================================================================##
italy_iucn_eval_edge <-  italy_iucn_eval %>% 
  mutate(CRB = ifelse(Category_CriteriaB == 'LC or NT', 
                      'NT', 
                      Category_CriteriaB)) %>% 
  filter(Taxon != 'Cycadeoidea_etrusca') %>% 
  dplyr::select(Taxon, CRB) %>% 
  rename(Redlists = CRB,
         Species = Taxon) %>% 
  as.data.frame() %>% 
  dplyr::select(Species, Redlists) 
##============================================================================##
italian_edge <- readRDS('Italian EDGE.rds')
##============================================================================##
## Estimate EDGE -----
##============================================================================##
italian_edge <- phyloregion::EDGE(italy_iucn_eval_edge, 
      italy_tree, 
     Redlist = "Redlists", 
     species = "Species")
##============================================================================##

##============================================================================##
## If the above throws an error, run the following chain of commands
##============================================================================##
# italian_edge <- phyloregion::evol_distinct(italy_tree,
#                                            type = 'fair.proportion') %>% 
#   enframe() %>% 
#   dplyr::rename(Taxon = name,
#                 ED = value) %>% 
#   inner_join(., italy_iucn_eval) %>% 
#   mutate(CRB = ifelse(Category_CriteriaB == 'CR', 0.999,
#                  ifelse(Category_CriteriaB == 'EN', 0.67,
#                         ifelse(Category_CriteriaB == 'VU', 0.1, 0.01))),
#          EDGE = EDGE_fnc(ED, CRB)) %>% 
#   dplyr::select(Taxon, ED, EDGE, everything())
# 
# map_edge <- italian_edge$EDGE
# names(map_edge) <- italian_edge$Taxon

##============================================================================##
##============================================================================##
## Create a duplicate for later analysis
##============================================================================##
italian_edge_tbl <- italian_edge %>% 
  enframe() %>% 
  dplyr::rename(Taxon = name,
                EDGE = value) %>% 
  inner_join(., clean_dts) %>% 
  dplyr::select(order, family, genus, Taxon, EDGE)
##============================================================================##
##============================================================================##
## Histogram
##============================================================================##
hist(italian_edge_tbl$EDGE, 
     main = 'Distribution of EDGE values in Italy', 
     xlab = 'EDGE')

##============================================================================##
## Basic statistical features
##============================================================================##
skimr::skim(italian_edge_tbl)


##============================================================================##
## Mean EDGE values by genus
##============================================================================##
italian_edge_tbl %>% 
  group_by(family, genus) %>% 
  summarise(Mean_EDGE = mean(EDGE),
            Species_Number = sum(NROW(genus))) %>% 
  arrange(desc(Species_Number)) %>% 
  ungroup() %>% 
  head()


##============================================================================##
## Find out which species have the top-1% of EDGE values
##============================================================================##
italian_edge_tbl %>% 
  arrange(desc(EDGE)) %>% 
  dplyr::select(Taxon, EDGE) %>% 
  filter(EDGE >= quantile(italian_edge_tbl$EDGE, 
         probs = 0.99, 
         na.rm = T))  %>% 
  head()
##============================================================================##
## Load the sparse matrix
##============================================================================##
pt <- readRDS('Sparse matrix and relevant polygons for Italy.rds')
##============================================================================##
##============================================================================##
## Create a sf object with hot-/cold-spots info
##============================================================================##
hot_cold_spots <- pt$poly_shp %>% 
  
  ## Convert to sf 
  st_as_sf() %>%
  
  ## Create the hot- and cold-spots
  mutate(Cold_SR = coldspots(richness, 
                             prob = 5),
         Hot_SR = hotspots(richness, 
                           prob = 5),
         
         ## Corrected weighted endemism
         Cold_CWE = coldspots(weighted_endemism(pt$comm_dat),
                              prob = 5),
         Hot_CWE = hotspots(weighted_endemism(pt$comm_dat), 
                            prob = 5),
         
         ## Phylogenetic diversity
         Cold_PD = coldspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                             prob = 5),
         Hot_PD = hotspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                              match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                           prob = 5),
         
         ## Phylogenetic endemism
         Cold_PE = coldspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                            match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                             prob = 5),
         Hot_PE = hotspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                          match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                           prob = 5),
         
         ## EDGE
         EDGE = map_trait(pt$comm_dat,
                          
                          ## If an error occured previously, change that with map_edge
                          italian_edge, 
                          FUN = median, 
                          shp = pt$poly_shp)$traits,
         Cold_EDGE = coldspots(EDGE, prob = 5),
         Hot_EDGE = hotspots(EDGE, prob = 5))
##============================================================================##
## Create a new variable
##============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% 
  as_Spatial()
##============================================================================##


##============================================================================##
## Quick plot
##============================================================================##
plot(sp_hot_cold_spots,
     border = "grey",
     col = "lightgrey",
     main = "Diversity Hotspots and Coldspots")

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Hot_EDGE == 1), ],
     col = "red",
     add = TRUE, 
     border = NA)

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Cold_EDGE == 1), ], 
     col = "blue", 
     add = TRUE, 
     border = NA)

legend("bottomleft",
       fill = c("red", "blue"),
       legend = c("hotspots", "coldspots"),
       bty = "n", 
       inset = .092)
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##

Βιβλιογραφία

Bachman, S.P., Nic Lughadha, E.M., & Rivers, M.C. (2018) Quantifying progress toward a conservation assessment for all plants. Conservation biology, 32, 516–524.

Daru, B.H., Roux, P.C. le, Gopalraj, J., Park, D.S., Holt, B.G., & Greve, M. (2019) Spatial overlaps between the global protected areas network and terrestrial hotspots of evolutionary diversity. Global Ecology and Biogeography, 28, 757–766.

Dauby, G., Stévart, T., Droissart, V., Cosiaux, A., Deblauwe, V., Simo-Droissart, M., Sosef, M.S., Lowry, P.P., Schatz, G.E., Gereau, R.E., & others (2017) ConR: An r package to assist large-scale multispecies preliminary conservation assessments using distribution data. Ecology and evolution, 7, 11292–11303.

Faith, D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological conservation, 61, 1–10.

Isaac, N.J., Turvey, S.T., Collen, B., Waterman, C., & Baillie, J.E. (2007) Mammals on the edge: Conservation priorities based on threat and phylogeny. PloS one, 2, e296.

Kougioumoutzis, K., Kokkoris, I.P., Panitsa, M., Trigas, P., Strid, A., & Dimopoulos, P. (2020a) Plant diversity patterns and conservation implications under climate-change scenarios in the mediterranean: The case of crete (aegean, greece). Diversity, 12, 270.

Kougioumoutzis, K., Kokkoris, I.P., Panitsa, M., Trigas, P., Strid, A., & Dimopoulos, P. (2020b) Spatial phylogenetics, biogeographical patterns and conservation implications of the endemic flora of crete (aegean, greece) under climate change scenarios. Biology, 9, 199.

Rosauer, D., Laffan, S.W., Crisp, M.D., Donnellan, S.C., & Cook, L.G. (2009) Phylogenetic endemism: A new approach for identifying geographical concentrations of evolutionary history. Molecular ecology, 18, 4061–4072.

Stévart, T., Dauby, G., Lowry, P., Blach-Overgaard, A., Droissart, V., Harris, D., Mackinder, B., Schatz, G., Sonké, B., Sosef, M., & others (2019) A third of the tropical african flora is potentially threatened with extinction. Science Advances, 5, eaax9444.

Tucker, C.M., Cadotte, M.W., Carvalho, S.B., Davies, T.J., Ferrier, S., Fritz, S.A., Grenyer, R., Helmus, M.R., Jin, L.S., Mooers, A.O., & others (2017) A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biological Reviews, 92, 698–715.


  1. Βασική προϋπόθεση φυσικά είναι να τις έχετε εγκαταστήσει ήδη

  2. Υπάρχουν δύο κατηγορίες ακόμα, όμως τα είδη τα οποία εμπίπτουν σε αυτές δεν αντιμετωπίζουν ακόμα εμφανή κίνδυνο

  3. Και μακροχρόνιαα εργασία πεδίου

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr